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ABSTRACT 



O 

o 

^ ■ We present a new "colhsional grooming" algorithm that enables us to model 

Q ■ images of debris disks where the collision time is less than the Poynting Robertson 

^ ■ time for the dominant grain size. Our algorithm uses the output of a collisionless 

disk simulation to iteratively solve the mass flux equation for the density distri- 
bution of a colhsional disk containing planets in 3 dimensions. The algorithm can 
|Vj ■ be run on a single processor in ~1 hour. Our preliminary models of disks with 

C/^ ■ resonant ring structures caused by terrestrial mass planets show that the colli- 

■ sion rate for background particles in a ring structure is enhanced by a factor of 

I . a few compared to the rest of the disk, and that dust grains in or near resonance 

2 ■ have even higher collision rates. We show how collisions can alter the morphol- 

^ . ogy of a resonant ring structure by reducing the sharpness of a resonant ring's 

inner edge and by smearing out azimuthal structure. We implement a simple 
prescription for particle fragmentation and show how Poynting- Robert son drag 
. and fragmentation sort particles by size, producing smaller dust grains at smaller 

CN . circumstellar distances. This mechanism could cause a disk to look different at 

^ . different wavelengths, and may explain the warm component of dust interior to 

Q>^ . Fomalhaut's outer dust ring seen in the resolved 24 fim Spitzer image of this 

O 

a^ 
o 



system. 

Subject headings: circumstellar matter — interplanetary medium — methods: 
N-body simulations — methods: numerical — planetary systems 



1. Introduction 

Recent resolved images of several debris disks reveal complex structures in the form of 



rings, gaps, and warps carved in the circumstellar dust (e.g. lGreaves et al.lll998l : IWilner et al. 
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20021 : iKalas et al.l l2005l : iGolimowski et all l2006l : ISchneider et al.l |2009| ) . Some of these struc- 
tures are likely the result of planeta ry companions that gravitationally perturb the disk 



[e.g. lQuillenll2006l : IChiang et al.ll2009l ). Modeling these structures can potentially reveal the 
physical and orbital parameters of the planets, dust grains, and sources of dust in these 



represent (e.g. 


Zuckerman 


2001; 


Wilner et al. 


2002; 


Moran et al. 


2004; 


Deller & Maddison 


2005; 


Wyatt et al. 


2007; 


Stark & Kuchnerl 


2008). 



So far all resolved debris disks are collisionally-dominated systems, meaning the c ollision 



time, tcoih is much shorter than the Poynting- Robertson (PR) time, tpR, in these disks ( IWyatt 
20051 ). For exa mple, tpR /tmn is ^3 00 for 10 /im grai ns in the resolved circumstellar dust ring s 



of Fomalhaut fIChiang et al.ll2009r ) and HR 4796A flDebes et al.ll2008l ; ISchneider et al.ll2009r ). 

Many different modeling techniques have been applied to these disks. However, no model 
has yet been able to accurately treat gravitational dynamics and collisions simultaneously 
in a self-consistent fashion. 



Some models of dust in debris disks ignore collisions altogether (e.g. iMoran et al.ll2004j ; 
Deller fc Maddison! |2005| ). Some variations on collisionless disk models simply stop the in- 
tegr ation of the particle orbits once dust grains have lived as long as their coUisional time 



[e.g. 



Chiang et al.ll2009l ). The collision time is typically estimated as tcou = ^orbit/aT, where 
^orbit is the orbital period, r is the optical depth of the disk, and a is a constant determined 
by the distribution of the grains' i nclinations and ecce ntricities and the size threshold for 
catastrophic collisions, estimated by lWyatt et al.l (Il999l ) to be a ~ in. These models do not 
include the influence of disk asymmetries and orbital resonances on collision rates. 

Other models of dusty disks avoid treating the resonant dynamics of the grains al- 
together, opting instead to investigate the long-term coUisional evolution of debris disks 
with simple geometries. Analytic models of azimuthally symmetric steady-state coUisional 



Wvatt 


1999; 


Wvatt et al. 


2007) 



with great detail by including processe s such as grain fragmentation and cratering (e.g. 



Thebault et al.ll2003l ; iLohne et al.ll2008l ). specialize in the long-term behavior of the radial 



an d grain size distribu t ions o f dust in disks, and do not include planets. The hybrid model 



of iBromely fc KenyonI (120061 ) merges kinetic models with n-body models by combining a 
multi-annulus kinetic model which treats small dust grains with an n-body model to treat 
large planetesimals. However, this code is tailored for grain growth simulations and not res- 
onant structures in debris disks; it does not model the resonant dynamics of the dust grain 
population. 

Still other models look at only the transient effects of collisions. The coUisional code 
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of iGrigorieva et al.l (120071 ) is designed for modeling collisional avalanche events in debris 
disks, but does not have the spatial resolution necessary to model structures caused by 
gravitational resonant dynamics. Robust many-particle simulations that model individual 
collisional events and follow the orbits of any fragments produced are curr ently limited by 



comp uter processing power to ~10^ particles for long-term integrations (e.g. iLeinhardt et al. 



20091 ). restricting their debris disk application to short-term integrations or integrations 
samphng limited phase space. 

Here we present a novel "collisional grooming" algorithm for treating collisions and 
resonant dynamics self-consistently in an optically thin disk. The algorithm produces a dust 
distribution for each grain size that simultaneously solves the equation of motion for small 
dust grains, 
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and the particle number flux equation in 3D, 
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Here G is the gravitational constant, is the stellar mass, c is the speed of light, r and 
V are the heliocentric position and velocity of the grain, /3 is the ratio of radiation pressure 
force to gravitational force on a grain, sw is the ratio of solar wind drag to PR drag, and rrii 
and Ti are the mass and heliocentric position of the z*^ planet. The particle number density, 
n, the mean flow velocity, v, and the particle number density removed or added by collisions 
per unit time, dn/dt\coii, niay be functions of position and grain size. 

Our algorithm uses a collisionless disk simulation (a "see d model" ) as input to c a lculat e 
initial collision rates, similar to the metho d developed bv ICharnoz fc Morbidellil (120031 ). 
However, our algorithm differs from that of ICharnoz fc Morbidellil (120031 ) in several ways. 
Our algorithm includes Poynting-Robertson drag and is therefore applicable to systems with 
small dust grains. Our algorithm also allows for the possibility that collisions can affect the 
dynamics of the system, and it uses an iterative scheme to find the correct density distribution 
for a steady-state disk. After the integration of the seed model, the algorithm can be run on 
a single processor in ~1 hour. This algorithm can generate new models that should allow us 
to interpret images of collision-dominated disks like those orbiting Fomalhaut, Vega, Epsilon 
Eridani, and HR 4796A quantitatively for the first time. 
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Numerical Method 



2.1. Collisionless seed models 



We first run a seed model, a model of a steady-state collisionless debris disk. We 
numerically integrate Equation [1] for a collection of particles launched from parent bodies 
orbiting a star. This equation includ es the dynamical effects of gravity, radiation pressure, 
corpuscular drag, and PR drag (see IStark fc Kuchneii l2008l . for details). The drag forces 
cause the particles to slowly lose angular momentum and spiral inward from the parent 
bodies from which they were launched. 

During their journey inward toward the star, the particles can become temporarily 
trapped in the external mean motion resonances (MMRs) of any planets present, creating 
an overdense circumstellar ring structure near the planet's orbit. We use enough particles, 
typically on the order of a few tho usand, to accurately populate the external MMRs of 
the planets (jstark fc Kuchne"r 2008 ). During the integration, we record the barycentric 
coordinates of each particle at a regular interval, ^record- This commonly- used technique 



extrapolates the results of a few thousand particles to millions of particles (e. g. iDermott et al. 
19941 : Ihiou k ZooklEoogl : iMoro-Martin k Malhotralbooi iMoran et al.ll2004h . Our algorithm 
requires the local velocity distribution to calculate the local collision rate, so we also record 
the barycentric velocities of each particle at the same interval, trecord- 

We then place all of the records of the barycentric coordinates and velocities into a 3D 
spatial grid of bins, forming a 3D histogram that represents the distribution function for the 
collisionless system. Panel b in Figure [1] shows the histogram of the particle density viewed 
face-on for a collisionless disk model with a resonant ring created by an Earth-mass planet 
on a circular orbit at 1 AU around the Sun, using a grid bin size of 0.05 AU. Panel a shows 
the velocity records from three of these bins located in the midplane of the disk, pointed to 
in panel b. 

Exterior to the resonant ring structure in Figure [H the velocity distribution is approx- 
imately Gaussian and the dispersion in the radial direction is approxir nately twice the dis- 
persi on in the azimuthal direction, as expected for a Keplerian disk (e.g. iBinney k Tremaine 
19871 ). However, within the resonant ring structure the velocity distribution is highly non- 
Gaussian and the velocity dispersion varies greatly from one location to another. This 
illustration shows that to calculate the collision rates in a resonant ring we must explicitly 
calculate the local velocity distribution. 
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2.2. The collisional grooming algorithm 

Besides the distribution function, recording the barycentric coordinates and velocities 
of all particles at regular intervals during the integration of the seed model yields a second 
important ingredient for our algorithm — a chronological record of each particle's trajectory, 
which we refer to as a "stream." We let each stream from our seed model represent a large 
number of particles, which varies from record to record, i.e. the i^^ record is scaled to Ni 
particles. We adjust the scaling factor to control the mean collision time of the particles. 

We initially assume collisions only serve to remove material from a stream. As they 
progress through the cloud, the streams become attenuated by collisions with other streams 
as 

iV, = iV,_i e-"-"'S (3) 

where Ni is the number of particles in the i^^ record of a given stream and rcoii,j is the 
collision depth for the i^^ record. Equation [3] is analogous to the solution to the radiative 
transfer equation for photons passing through an absorptive medium. We approximate the 
collisional depth as 

T"coll,i ~ ^^fcO-fclVj - Vfc|trecord, (4) 

k 

where Vj is the velocity associated with the i^^ record of the given stream, n^, ak, and 
are the particle number density, collisional cross-section, and velocity of the other records 
in the same bin as the i^^ record, and is equal to divided by the bin volume. This 
approximation works as long as trecord tcoii- We perform this calculation for all records in 
all particle streams, one at a time, from the first record to the last record in each stream. 

After only one pass through all of the records for all streams, i.e., one iteration, the 
particle streams will be attenuated incorrectly; the streams will be attenuated based on the 
density distribution of the coUisionless seed model, which overestimates the particle density 
and therefore the collision rate. To remedy this problem, we iterate the attenuation process 
of all records until no record changes by more than a set tolerance of a few percent. By doing 
so, we ensure that the final number density histogram approximately satisfies Equation [2] in 
the steady state. 

Figure [2] shows an example of the grooming algorithm at work. The top-left panel shows 
the surface density of a coUisionless disk scaled by where r is circumstellar distance. 

The disk of 120 fim grains features a ring structure caused by an Earth-mass planet on a 
circular orbit at 1 AU around a Sun-like star. We applied our collisional grooming algorithm 
to the disk and scaled the number of particles per stream equally among all streams such 
that ?7o ~ 3.7, where rjo = tpR(ro)/tcoii('"o) and tq is the mean circumstellar distance at which 
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grains are launched. 

After the first iteration, shown in the top-middle panel, the collision rates are overesti- 
mated so that the surface density in the inner disk is too low. During the second iteration, 
shown in the top-right, the algorithm underestimates the collision rates. The algorithm alter- 
nates between over- and underestimating the collision rates while converging on the correct 
solution. 

Our algorithm uses a finite grid of bins to approximate the local density and velocity 
structure, so Poisson noise can become an issue in bins with few records. To help mitigate 
this Poisson noise, we use a nearest-neighbor averaging routine. For bins with fewer than 10 
records, we include the records from the six nearest-neighbor bins in the calculation of Tcou. 
We weight the records in each nearest-neighbor bin at 50% of the weight of the central bin. 

Our calculation of the collisional depth (Equation Hj) trivially handles seed models with 
multiple grain sizes. Given an initial distribution of grain sizes, the algorithm can simulta- 
neously solve for the collisional interactions among all grains of all sizes included in the seed 
model. The algorithm can model complex phenomena that depend on grain size, such as 
size-dependent collision rates, radial transport rates, and resonant structure morphologies. 

Equation [3] ignores any fragments that may be produced by collisions. Treating collisions 
as non-productive is not valid for all collisions in all debris disks, but it is likely acceptable 



in a wide range of cases. iBackman &: Parescd (119931 ) argued that in the inner solar system 



grain-grain collision velocities are high enough (on the order of a few kilometers per second) 
that grains fragment catastrophically during any collision, and any fragmentation products 
are either gas or small enough to be removed immediately by radiation pressure. In resonant 
ring structures, like the circumsolar ring created by Earth, the collision velocities can be even 
higher, on the order of tens of kilometers per second (c.f. Figured]). In the outer solar system, 
where collision velocities are lower, grains probably resemble cometary grains. Samples 
of cometary interplanetary dust particles directly returned from the Stardust mission and 
observations of cometary ejecta during the Deep Impact mission reveal that the majority of 
observed cometary particles are loosely bound aggreg ates of submicron-siz e d grains, which 



can easily be shattered in to unbound /5-meteoroids (lA'Hearn et al.l l2005l : iBrownlee et al. 



200a IZolenskv et al.ll2006h . 



Our algorithm can also be adapted to handle particle fragmentation. We discuss this 
feature in Section HI 
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2.3. Tests 

We developed a collisional disk modeling code based upon the algorithm described 
above. We subjected our code to a battery of tests to confirm its operation and identify 
its limits. Here we demonstrate the algorithm's performance, show that it converges on a 
unique and correct solution, and place limits on the conditions under which it is valid. For 
now, we neglect fragmentation; we assume non-productive collisions. 

Below, we will refer to an analytic sol ution f o r the surface density as a function of 



circumstellar distance for a planet-less disk. IWyattI (119991 ) showed that Equation [2] has the 



following steady-state solution for an azimuthally-symmetric disk with a single grain size 
and grain mass under the assumption that collisions create no daughter particles: 

S(r) = ^ , (5) 

1 + 4r/o(l - Vr/ro) 



where E is the surface density and ro is the circumstellar distance of the dust source. If we 
do not add a perturbing planet, we can directly compare the results of our algorithm to this 
expression. 



2.3.1. Seed models 

Throughout this paper, we will refer to two collisionless disk simulations which we use 
as seed models: a planet-less disk simulation and a simulation of a disk with a resonant 
ring structure. For both simulations, we integrated the orbits of 20,000 particles. Particle 
integrations were terminated when their semi-major axes were less than 0.5 AU, so many 
images will show no data interior to this radius. In many instances throughout this paper 
we will only examine a subset of the total number of integrated particles, and will state the 
number used when appropriate. All s imulations were performed using a hybrid symplectic 



integrator (see lStark fc Kuchnerl |2008| ) . 



Our planet-less disk seed model contains particles with P = 0.0023, where P is the ratio 



of th e force on the grain from radiation pressure to the gravitational force (e.g. iBurns et al. 
19791 ). For a spherical blackbody grain with a density of 2 gm cm~^, (3 = 0.0023 corresponds 



to a grain radius of 120 /im around a Sun-like star. We initially placed the grains at 10 AU 
on circular orbits with inclinations uniformly distributed between and 14°. We distributed 
all other initial orbital parameters (longitude of ascending node, mean anomaly, argument of 
pericenter) uniformly between and 271. We recorded particle positions and velocities every 
6956 years. 
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For the second seed model, we integrated the orbits of particles with (3 = 0.0023 as 
they spiraled inward and interacted with an Earth-mass planet on a circular orbit at 1 
AU around a Sun-like star. We launched grains from parent bodies with semi-major axes 
uniformly distributed between 3.5 and 4.5 AU, eccentricities uniformly distributed between 
and 0.2, and inclinations uniformly distributed between and 20°. We distributed all other 
initial orbital parameters uniformly between and 27r. We recorded the particle positions 
and velocities once every 426 years. 



2.3.2. Bin size test 

The finite size of the bins used to approximate the local density and velocity distributions 
is a natural source of error for our algorithm. We need to ensure that the bins are small 
enough to resolve any structure within the disk. To help us decide on the appropriate bin 
size, we performed the following test. 

We applied our collisional grooming algorithm to our seed model of a disk with a resonant 
ring structure. We used 15,000 simulated particles to ensure that Poisson noise was not an 
issue for this test. For the collisional algorithm, we scaled the number of particles per stream 
such that riQ ~ 1, and used four different cubic bin sizes of 0.02, 0.05, 0.1, and 0.2 AU. For 
each of the largest three cases, we calculated the differences in the collision rates compared 
to the smallest case. 

Figure [3] shows these relative differences in the collision rates. The left panel shows that 
differences in the collision rates between the 0.05 and 0.02 AU bin sizes are on the order 
of a few percent or less, and show no signs of structural differences. The right panel shows 
that using a bin size of 0.2 AU results in obvious structural differences (greater than 10%) 
compared to a model with a 0.02 AU bin size. The right panel also shows a subtle imprint 
of the grid itself, which appears as straight horizontal and vertical features in the collision 
rates. 

In light of these results, we recommend that the bin size for collisional calculations in 
a debris disk with a resonant ring structure should be ~ 0.05 AU for a planet at 1 AU. 
The size of structural featu res in a resonant ring scale linearly with planet semi-major axis. 



Op (IStark fc Kuchnen 120081 ) . so we suggest that in general, the optimal bin size for disks 
with resonant ring structures is ~ 0.05 Op. Bins larger than this size fail to resolve the ring 
structure while bins smaller than this size become more susceptible to Poisson noise in the 
distribution function. 



For collision times that are short compared to the PR time [tiq ^ 1), grains launched 
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from circumstellar distance tq have little time to move radially inward before colliding. This 
process can result in a very narrow ring structure at r = Tq. If the bin size is larger than 
the width of the ring, our coUisional algorithm will fail to resolve the ring. So we must also 
choose a bin size small enough such that the crossing time of a bin caused by PR drag at 
r = ro is shorter than the collision time. The bin crossing time is given by 

Across (ro) = tpR(ro) - tpR(ro - h) , (6) 
where b is the bin size and the PR time is given by 

2 

where c is the speed of light, G is the gravitational constant, and is the stellar mass 



(IWyatt &: Whipplelll95Cll ). With the requirement tcross('^o) ^ ^coul'^o) we find to first order in 



^ that the bin size must satisfy 



2.3.3. Poisson noise test 

If there are too few records in a given bin, Poisson noise can dominate the calculation of 
the collisional depth, even with our nearest-neighbor averaging routine. However, we are not 
specifically concerned with whether a single collisional depth calculation suffers from Poisson 
noise, but whether Poisson noise affects the final outcome of the simulation. To examine the 
effects of Poisson noise, we applied our collisional algorithm to our planet-less seed model 
three times, first using 15,000 simulated particles, then using 5,000 simulated particles, and 
then using only 1,500 simulated particles. For each application of the collisional algorithm, 
we used a cubic bin size of 0.05 AU. We scaled the disk such that r^o ~ 223 for all three 
cases. 

After processing the three disks with our collisional grooming algorithm, we azimuthally 
averaged each of the three resulting surface densities. Figure H] shows the normalized surface 
density as a function of circumstellar distance for all three cases compared to the analytic 
solution (Equation [5]). The 15,000-particle case follows the analytic solution well, while the 
1,500-particle case deviates by a factor of ~ 2 and contains strong fluctuations caused by 
Poisson noise. From these simulations we estimate that the average number of records per 
bin near the circumstellar distance at which grains are launched should be at least on the 
order of a few to avoid the effects of Poisson noise. 
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2.3.4- Uniqueness of solution 

Figure [2] shows that our algorithm does indeed converge, but does not indicate whether 
the algorithm converges on a unique solution. To test for uniqueness, we applied the coUi- 
sional grooming algorithm to three independent seed simulations of disks with ring structures, 
using 5,000 particles for each simulation and a value of r/o ~ 4. The range of initial condi- 
tions were the same for all three simulations, but individual values were generated using a 
different random number seed. 

We compared the resulting surface densities of all three collisional disks. Except for 
the outer and inner extremities of the disk where the statistics were poor, none of the final 
disk surface densities differ by more than a few percent and none show significant structural 
differences. All three independent simulations converged to the same surface density solution 
to within the limits of Poisson noise. 



2.3.5. Correctness of solution 

To test the correctness of our algorithm's solution, we directly compared the results of 
our algorithm to the analytic solution for a planet-less disk of a single grain size (Equation [5]). 
We applied our collisional grooming algorithm to our planet- less-disk seed model of 20,000 
particles using six different values of r/o- Figure [5] shows the azimuthally-averaged surface 
densities as a function of circumstellar distance for all six cases. The calculated surface 
densities (solid lines) match the analytic solutions (dashed lines) well for all values of r/o 
except rjQ = 763. In this case, collisions happen so quickly that the disk forms a narrow ring 
at 10 AU whose width is smaller than the bin size used in our algorithm (0.05 AU); this case 
does not meet the bin size criteria in Equation [81 If we used a smaller bin size, we could 
resolve the ring structure and investigate even larger values of tjq. 

For all cases shown in Figure O trecord ^ ^cou; none of the deviations from the analytic 
values are caused by insufficient time sampling. We tested the algorithm's behavior when 
^record ^ ^coib a^d Confirmed that our algorithm fails under these circumstances. For trecord ^ 
tcoib our algorithm typically overestimates the collision rate and the result looks qualitatively 
similar to the tjq = 763 case shown in Figure [51 

For the case of a disk with a resonant ring structure, there exists no analytic solution 
for the surface density with which we can compare the results of our simulations. However, 
we can probably assume that our collisional algorithm arrives at the correct solution if the 
amount of collisions is very small, i.e., rjQ <^ 1, since such a disk should deviate from the 
collisionless case by very little. Under this assumption, we propose the following test to 
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investigate the correctness of our algorithm's solution for a disk with a ring structure: 

1. Apply our algorithm to the collisionless seed model using rjo 1. 

2. Store the output of the algorithm, call it Model C. 

3. Increase rjo by a small amount 6riQ. 

4. Apply our algorithm to Model C disk using the new value of rjo. 

5. Repeat steps 2 through 4 until r^o is equal to the desired value of tiq. 

6. Compare the results to a model with the same value of rjQ calculated in only one step. 

We performed this test using a 5,000-particle seed model of a disk with a ring structure. 
We used a cubic bin size of 0.05 AU and scaled the disk density so that the final rjo ^ 3.7. 
We applied the coUisional grooming algorithm using 19 logarithmically-spaced steps in tjq. 
We compared the surface density of this disk to the surface density calculated by applying 
the algorithm in the usual single-step fashion. The disks differed by less than one part in 
~ 10^; both methods arrived at the same solution to within the limits of Poisson noise. A 
similar test performed in the opposite direction, i.e., slowly reducing the rjo to the desired 
value, gave similar results. This test supports both the uniqueness and correctness of the 
solution that our algorithm finds for a coUisional disk with a resonant ring structure. 

2.3.6. Benchmark tests 

We benchmarked our coUisional grooming algorithm code by applying it to our seed 
model of a disk with a ring structure. We recorded the run time of our code using 1,250, 
2,500, 3,750, and 5,000 simulated particles and cubic bin sizes of 0.05 and 0.1 AU. Tabled] 
shows the run time for a single iteration of our algorithm on a single 2.2 GHz CPU. 

For each bin, the algorithm performs of order calculations, where nb is the number 
of records in that bin. So we would expect that our run time per iteration scales as B{n\,)'^, 
where B is the number of bins containing records, or as (^b)^ for a given bin size. Tabled] 
shows that the run time scales as (nb)^ for a bin size of 0.1 AU, but not for a bin size of 
0.05 AU. In the latter case, our algorithm is working with many bins that have relatively 
few entries and is switching on our nearest-neighbor approximation described in Section 12. 2[ 
which can cause the run time to deviate from the (^^b)^ scaling relationship. 

The number of iterations required to converge on the correct solution depends on many 
factors, such as r^o; the number of records, the bin size, the structure of the disk, and the 
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tolerance set for convergence. For a disk with a resonant ring structure composed of 5,000 
particles, typically 5-10 iterations are required for better than 5% convergence with r/o ~ 1- 
The total run time for such a disk tj^ically ranges from 20 minutes to 2 hours. 



2.4. Limitations 



Here we summarize the current limitations of our algorithm. We have already shown 
that the bin size must be small enough to resolve any structural features in the disk, and the 
collision time must be longer than the time between records for our algorithm to converge. 
We also showed that the number of records per bin must be large enough to avoid the effects 
of Poisson noise. 

The algorithm presently has another, more s ubtle limitation, which it shares with 



Dermott et al. 


1994; 


Liou & Zook 


1999 



Moran et al.ll2004l : iDeller fc Maddisorul2005l ). As stated in Section [221 we record the coordi- 
nates and velocities of each integrated particle at regular time intervals, trecord, to extrapolate 
the results of a few thousand integrated grains to millions of grains or more and to obtain a 
chronological record of each particle's trajectory. This technique implies that a new parent 
body launches grains every trecord with the exact same orbital parameters as the first parent 
body. Because our algorithm requires the use of this technique, our algorithm implicitly 
assumes that each parent body's orbit is populated by many parent bodies uniformly dis- 
tributed in mean anomaly and that the orbits of parent bodies do not change over the course 
of a coUisional time. 



3. Non-productive collisions in resonant ring structures 

We can use our algorithm to investigate the effects of collisions on steady-state resonant 
ring structures. To this end, we ran some models of familiar kinds of disk structures. Here 
we discuss some of the new physical effects we have observed in our models. 

The top-left panel of Figure E] shows the collision rate per particle in a 0.4 AU-thick 
cross-section through the mid-plane of a disk. The disk has a resonant ring structure caused 
by an Earth-mass planet whose location is marked with a white dot. We scaled the disk 
density so that tiq ^ 3.7. The top-right panel shows the surface density of the same cross- 
section for comparison. 

The collision rate is high in a circumstellar ring near the location of the parent bodies 
(~3 AU), i.e., the birth ring. Grains in this region of the disk are too young to have been 
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destroyed by collisions, so the local density is relatively high, as seen in the top-right panel. 
The collision rate drops just interior to this ring, at a circumstellar distance of ~2 AU. 

As you might expect, the top-left panel of Figure [H] shows that the collision rate reaches 
its highest point in the resonant ring structure, where the density is enhanced by resonant 
trapping and relative velocities are higher because of resonant pumping of the grains' ec- 
centricities. The average collision rate in the resonant ring structure is higher by a factor 
of a few when compared to the collision rate exte rior to the resonant ring structure, con- 



sistent with analytic estimates (iQueck et al.l 120071 ). Figure |6] also shows that the collision 



rate exhibits azimuthal and radial structure in the resonant ring. This structure reflects 
both localized density enhancements and regions of higher velocity dispersion. For example, 
the region of enhanced collision rate located ~90° clockwise from the planet generally corre- 
sponds to a region of higher density. But the region of enhanced colhsion rate located ~90° 
counterclockwise from the planet does not. This second region of enhanced collision rate is 
primarily caused by an increase in the local velocity dispersion, as shown in Figure [H 

We show the collision rate and surface density for an edge-on cross-section of the same 
disk in the bottom two panels of Figure El The bottom-left panel reveals a trend toward 
higher collision rates in the disk mid-plane, which is denser than the rest of the disk, as 
shown in the bottom-right panel. The bottom-left panel of Figure [6] also shows that the 
collision rate at a circumstellar distance of ~0.7 AU is higher than the collision rate at 
~2 AU by a factor of ~2, even though the density is higher near ~2 AU. This increase 
in the collision rate occurs because grains that survive and spiral inward past the resonant 
ring structure have typically had their eccentricities pumped up by passage through the 
resonance, so the velocity dispersion interior to the resonant ring structure is higher than 
the velocity dispersion exterior to the resonant ring structure. 

The top panel in Figure [7] shows the collision rate per particle as a function of semi-major 
axis in the region of the resonant ring, a kind of continuum with spikes. The continuum of 
the collision rate is higher for semi-major axis values close to 1 AU than for larger semi-major 
axis values because the resonant ring enhances the local collision rate by a factor of ~ 2, as 
shown in Figured The spikes in the collision rate correspond to MMRs: the 2:1, 5:3, 3:2, 
7:5, 4:3, 9:7, 5:4, etc., from right to left. Most of the first order resonances (p+l:p) show a 
split peak, with higher collision rates at the edges of the resonance than at the center. The 
split peaks may be caused by collisions between the grains just outside of the MMRs with 
the grains in the MMRs. 

The dashed line in the top panel of Figure [7] shows the classically-calculated collision 
rate, l/tcou = otT /torhit- We set a = 2.6 so that the collision rates w ere equal at a semi-m ajor 



axis of 1.5 AU. This value of a is less than the estimate a ~ 47r by lWyatt et al.l (Il999l ). and 
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hence our co l lision rate is lower, likely because we simulate only a single grain size whereas 



Wyatt et al.l (Il999l ) consider a distribution of grain sizes and include catastrophic collisions 



caused by projectile grains much smaller than the target grain. 

Our algorithm shows how the classically-calculated collision rate fails for disks with 
planets, since this approximation neglects collisional enhancement of particles in resonance 
and does not accurately approximate the average collision rate in the resonant ring structure. 
It also cannot correctly reproduce the vertical structure of the collision rate, like that shown 
in the bottom- left panel of Figure O The failure of this approximation is even more dramatic 
than these figures show because the particles spend most of their time in the MMRs. 

The bottom panel of Figure [7] shows the collision rate per particle as a function of 
eccentricity. The particles attain an eccentricity of no more than e ~ 0.2 when they are 
launched. So any particles with e > 0.2 in this plot must have had their eccentricities 
increased by resonant pumping or close encounters with the planet. These particles have 
higher collision rates because they are typically located in the resonant ring, a region of higher 
density and larger velocity dispersion. The small peaks in the collision rate at e = 0.27, 
e = 0.31, and e = 0.37 correspond to t he maximum eccentricities th at particles can obtain in 



the 5:4, 4:3, and 3:2 exterior MMRs (IBeauge fc Ferraz-Melldll994j ). and also correspond to 



localized regions of higher density within the resonant ring structure. The data become noisy 
for e > 0.38 because relatively few particles end up with such large eccentricities, except in 
close encounters with the planet. 

Figure [8] illustrates some of the morphological effects collisions can have on resonant 
ring structures. The top row shows surface density histograms of a disk with a resonant ring 
structure for three different values of rjo, increasing from left to right. The bottom row shows 
zoomed-in views of the resonant ring structures from each of the disks in the top row. Each 
of the six histogram images was scaled independently to highlight differences in geometry. 

Collisions de-emphasize the resonant ring and emphasize the birth ring. They also 
change the morphology of the resonant ring. In the coUisionless case on the left, the resonant 
ring structure has a sharp inner edge with well-defined azimuthal and radial structure. For 
the slightly collisional system [tiq ~ 0.09), shown in the middle, the density at the resonant 
ring's inner edge is significantly reduced relative to the rest of the ring structure. As we 
further increase the overall collision rate in the disk, the density of the inner edge of the 
resonant ring continues to drop relative to the rest of the ring, while azimuthal structures in 
the ring become smeared out. 



Collisions reduce the density of the inner edge of the resonant ring for two reasons. 
First, as previously shown, collision rates are higher in regions of higher density. Second, 
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grains that contribute to the inner edge of the resonant ring are typically older, having had 
their eccentricities pumped up while in resonance. These older grains have had more time 
to collide with other grains. 



4. Particle fragmentation 

Our algorithm, as described above, can handle seed models with multiple grain sizes with 
no modifications. With a small modification the algorithm can also model fragmentation, 
i.e., collisions that produce daughter particles. Here we describe a method for including 
particle fragmentation and present some preliminary results. 

When two particles collide and produce fragments, the fragments are launched into 
new orbits such that the fragment velocity vectors are distributed around the center of 



momentum velocity vector of the colliding particles (e.g. iKrivov et al.ll2005l ). Integrating 
the orbits of all of these fragments would be computationally prohibitive, so we make a 
numerical approximation: we limit the trajectories of the fragments we model to the recorded 
trajectories that already exist in the seed model. To make this approximation work, our seed 
model must include a sufficiently wide range of initial conditions to ensure that trajectories 
are available that closely match the desired distribution of fragment trajectories. 

To include particle fragmentation in our algorithm, we first run a seed model with several 
discrete particle sizes, each of which represents a range o f sizes, and initi ally populate each 



size bin according to a mass distribution function (e.g. iDohnanyil Il969l ). Then we apply 
our collision algorithm with the following additional subroutine, implemented during the 
calculation of the collision depth (Equation H]) for every record of each stream: 

1. Calculate the k^^ record's contribution to the collision depth of the i**^ record (c.f. 
Equation H]), call it 

Tcoll,i,A; = nkakl^i — Vfe|trecord- (9) 

Remember that the index i refers to records in a stream while the index k refers to 
records in a bin. 

2. Calculate the mass of particles removed from the i^^ record by collisions with the k^^ 
record, given by 

AMi,fc = m,Ar,_i (1 - e-^-"-^'=) , (10) 
where rrii is the mass of a single particle in the i^^ record. 



3. 



Record the center of momentum velocity vector of the colhding particles, given by 



VC0M,j,A: — 



rrii + ruk 



(11) 



We assume the difference between the fragment velocities and vcoM,j,fc is small and use 
vcoM,i,fc as the desired fragment velocity. 

4. Distribute the fragments by size according to a crushing law. 

5. Search within the local spatial bin for streams that closely match the grain size, s, and 
center of momentum velocity vector, vcoM,i,fc, of the fragments. Once the appropriate 
streams are found, increase the numbers of particles in those streams to account for 
the fragments. 

The subroutine described above must be executed for every pair of records in every 
bin. In practice, searching for an appropriate stream in which to put fragments during every 
interaction is computationally expensive; the algorithm run time scales as (nb)^ instead of 
(rib)^. One could imagine many possible approximations that would reduce the amount of 
computer time spent searching for fragment streams. For our preliminary models, we chose 
to place all of the fragments from the i^^ record into a single mean center of momentum 
stream, with velocity 



k 

We leave a detailed investigation into the accuracy and efficiency of this approximation for 
future work. 

We have implemented this particle fragmentation subroutine in our code and produced a 
simple model of a fragmenting disk with a resonant ring structure to illustrate the procedure 
at work. For our seed model, we integrated the orbits of 2,500 120 /xm grains (/5 = 0.0023) 
and 2,500 12 /im grains (/5 = 0.023) and recorded the particle positions and velocities every 
426 and 42 years, respectively, as they orbited a Sun-like star (sw= 0.35) in the presence 
of an Earth-mass planet on a circular orbit at 1 AU. We launched the grains from parent 
bodies with initial conditions identical to the second seed model described in Section 12.3.11 
For the purposes of illustration, we initially populated only the 120 /xm grain streams and 
implemented a simple fragmentation scenario in which 120 /im grains shatter completely 
into 12 ixm grains, conserving mass. Fragments from 12 /xm grains were deleted from the 
model. We assumed all colliding grains were shattered completely, regardless of the mass or 





k 



(12) 
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velocity of the target or projectile. We scaled the number of 120 /im particles per stream 
such that riQ ~ 0.04. 

Figure M shows the results of our fragmenting disk model. The inset false-color image 
shows the face-on surface density of the disk. The 120 fim (large) grains, shown in red, 
dominate the surface density exterior to the resonant ring structure. The 12 fim (small) 
grains, shown in blue, dominate the surface density near the center of the disk. 

The plot in Figure M shows the disk mass distribution as a function of semi-major axis 
for each of the grain sizes in our simple model. At the outer edge of the disk, near the 
birth ring of large grains at ~4 AU, the large grains dominate the mass of the disk. As the 
large grains spiral inward, coUisional fragmentation reduces the disk mass in large grains and 
transfers that mass to the small grains. 

Near 1.5 AU, the resonant ring structure enhances the collision rates of the large grains 
and therefore also the disk mass in small grains. The spikes in the mass distribution function 
near 1 AU show that the two grain sizes populate different sets of MMRs. The combined 
effects of collisional fragmentation and PR drag would cause the resonant ring structure, and 
the disk as a whole, to look different at different wavelengths, because grains of different sizes 
emit differently. Figure [9] shows that PR drag and resonant interactions can sort collision- 
ally fragmenting grains by size, allowing smaller grains to spiral in to smaller circumstellar 
distances. 

The radial distribution of dust grains produced b y our model is analogous to that 



observed in the Fomalhaut disk. IStapelfeldt et al.l (120041 ) resolved the Fomalhaut disk at 24 



70, and 160 /im using the Multiband Imaging Photometer for Spitzer (MIPS) and obtained 
a 17.5-34 /im spectrum with the Infrared Spectrograph (IRS). Both the IRS spectrum and 
the 24 yum MIPS image reveal a compact source of dust likely interior to 20 AU, well inside 
of Fomalhaut 's outer ring structure near 140 AU. This compact source of dust, responsible 
for ~0.7 Jy of flux at 24 /im, is not seen in the 70 or 160 /im MIPS images, suggesting that 
the warm dust grains are inefficient at radiating at these wavelengths; the compact warm 
dust grains may be smaller in size than the grains near the outer ring. Our preliminary 
fragmentation model shown in Figure suggests that collisional fragmentation of large grains 
triggered by MMRs may be the source of small dust grains. 



5. Summary 



We have developed a new algorithm to self-consistently treat collisions and resonant 
gravitational dynamics in dusty disks. Our algorithm handles disks with multiple grain sizes 
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and can be adapted to model particle fragmentation. The algorithm uses the density and 
velocity distributions of a collisionless disk simulation to iteratively solve for the density 
distribution of a steady-state collisional disk. The algorithm is applied after the simulation 
of the collisionless system, removing the need to re-integrate the equations of motion for 
disks with different collision rates, and can run on a single processor in ~ 1 hour. 

We performed several tests to show that our algorithm arrives at a unique and correct 
solution for collisional disks with and without a resonant ring structure. We showed that 
collisions can reduce the contrast of resonant ring structures, especially at the inner edge 
of the ring structure, and smear out azimuthal asymmetries. We also showed that particle 
fragmentation triggered by resonant interactions can radially sort particles by size, producing 
smaller particles at smaller circumstellar distance s. This process may exp lain the population 



of warm dust found interior to Fomalhaut's ring (jStapelfeldt et al.ll2004l ) 



Our collisional grooming algorithm should allow us to accurately model and synthesize 
multi-wavelength images of observed debris disks, like Fomalhaut, Vega, and HR 4796A. The 
algorithm enables us to investigate the effects that collisions have on dust disk morphology, 
such as asymmetries from clumping of parent bodies, resonant trapping of dust grains, and 
the radial sorting of grain sizes illustrated in Figure [91 It should be useful for modeling 
long-lived grains in the solar zodiacal cloud and it should help us predict the morphology of 
ring structures in disks yet to be observed. 
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Fig. 1. — Velocity distributions (a) at three locations (6) in a collisionless dust disk with a 
resonant ring structure caused by an Earth-mass planet at 1 AU. The velocity distributions 
in the resonant ring structure vary greatly and are non-Gaussian. 
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Fig. 2. — Surface density scaled by of a collisional model as a function of iteration 
number. The upper-left panel shows the surface density of the seed model, a coUisionless 
disk with a resonant ring structure (see Section [2^3] for simulation details). The iterations 
alternate between over- and underestimating the amount of collisions while converging on 
the correct solution. 
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Fig. 3. — Errors in the collision rate in a model with an Earth-mass planet at 1 AU. Errors 
are relative to the collision rates calculated using a bin size of 0.02 AU. 
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Fig. 4. — Azimuthally averaged surface density as a function of circumstellar distance in 
a planet-less collisional disk for three simulations with different numbers of particles (see 
Section [2^3]) ■ The 15,000-particle simulation follows the analytic solution well, but the 1,500- 
particle simulation deviates significantly and shows signs of Poisson noise. 
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Fig. 5. — Azimuthally averaged steady-state surface density as a function of circumstellar 
distance for a planet-less disk undergoing PR drag and non-productive collisions. Analytic 
solutions are shown with dashed lines and the results of our collisional algorithm are shown 
with solid lines. Our algorithm gives the correct solution for all values of r]o except the largest, 
at which point the collisions yield a narrow ring near 10 AU that is no longer resolved by 
the grid. 
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Tablc 1. Benchmark Tests 



Number of Particles 


Bin Size 


Run Time per Iteration^ 
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0.05 
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0.05 
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0.05 
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1,250 


0.1 
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0.1 


4.2 
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0.1 


16.7 



tPor a single 2.2 GHz CPU 
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Fig. 6. — Top-left: Collision rate per particle {\dN/Ndt\) in a 0.4 AU-thick cross-section 
through the midplane of a disk with a resonant ring structure. A white dot marks the 
location of an Earth-mass planet orbiting at 1 AU. Top-right: Surface density of the midplane 
cross-section shown in the upper- left panel. Bottom-left: Collision rate per particle in a 0.4 
AU-thick edge-on cross-section of the same disk. Bottom-right: Surface density of the edge- 
on cross-section shown in the bottom-left panel. The collision rate is affected by both the 
local density structure and the local velocity distribution; the collision rate is highest in the 
resonant ring structure, a region of enhanced density and relative particle velocities. 
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Fig. 7. — Collision rate per particle {\dN/Ndt\) as a function of semi-major axis and eccen- 
tricity for a disk with a resonant ring structure. The spikes in the collision rate vs semi-major 
axis and eccentricity show that the collision rate is enhanced for particles in resonance, and 
also adjacent to the first order MMRs. 
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Fig. 8. — Surface density as a function of rjo for a disk with a ring structure caused by an 
Earth-mass planet at 1 AU viewed face-on. The top row shows the entire disk, which extends 
out to 4.25 AU. The bottom row shows zoomed-in views of the resonant ring structure. 
CoUisions reduce the sharp inner-edge feature of the resonant ring structure, smear out 
azimuthal structure, and de-emphasize the resonant ring while emphasizing the birth ring. 
Even a low collision rate {t]q '^1) can significantly alter resonant ring structures. 
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Fig. 9. — Disk mass (in Lunar masses) as a function of semi-major axis for a disk of frag- 
menting grains in the presence of an Earth-mass planet at 1 AU orbiting the Sun. The inset 
false-color image shows the face-on surface density of the disk. MMRs near 1 AU trigger 
fragmentation, a process which may explain the population of small warm dust interior to 
Fomalhaut's resolved ring structure (jStapelfeldt et al.ll2004j ). 



